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ABSTRACT 


A  brief  summary  of  the  research  efforts  under  this 
contract  is  provided.  The  areas  of  investigation  fall  under 
five  categories:  1)  Seismic  wave  propagation  in  hetero¬ 
geneous  media  where  finite  difference  and  other  numerical 
techniques  are  applied  to  non-uniform  geometries,  2)  Seismic 
source  theory,  and  the  determination  of  source  parameters 
and  focal  depth  from  spectra  of  Rayleigh  and  Love  waves, 

3)  Study  of  structure,  properties  and  heterogeneities  of 
the  earth's  interior  from  seismic  body  and  surface  waves  and 
generation  of  average  earth  models,  4)  Utilization  of 
seismic  arrays,  and  automatic  schemes  of  phase  identification 
and  event  verification  from  the  LASA  detection  logs, 

5)  Development  and  performance  of  long  period  tiltmeters. 

A  comprehensive  list  of  articles  published  on  the  above 
subjects  is  provided. 


11 


TABLE  OF  CONTENTS 

I.  INTRODUCTION 

II.  WAVE  PROPAGATION  STUDIES 

III.  SOURCE  MECHANISM  AND  FOCAL  DEPTH  STUDIES 

IV.  EARTH  STRUCTURE  AND  PROPERTIES 

V.  SEISMIC  ARRAY  STUDIES 

VI.  INSTRUMENTATION 

List  of  Publications 
List  of  Theses 


PAGE 

1 

2 

7 

11 

19 

49 

54 

60 


1. 


I.  INTRODUCTION 

In  this  final  report  we  will  briefly  summarize  the 
work  done  under  the  contract  during  the  past  five  years. 

As  the  title  of  the  contract  indicates,  this  project  was 
intended  as  a  study  of  broad  aspects  of  seismology  and  the 
understanding  of  the  specific  problems  within  the  general 
development.  The  studies  cover  the  elastic  and  non-elastic 
properties  of  the  earth  and  the  nature  of  earthquake  and 
explosive  sources.  Greatest  efforts  were  made  to  go  beyond 
the  classical  techniques  in  studying  general  problems. 

The  areas  of  investigations  can  be  combined  under  five 
general  headings: 

1.  Seismic  wave  propagation  studies 

2.  Source  mechanism  and  focal  depth 

3.  Earth  structure  and  properties 

4.  Seismic  array  studies 

5.  Instrument  developments 

Because  of  length  limitations,  we  cannot  discuss  in  detail 
all  the  work  done  under  each  heading.  Instead  we  introduce 
a  brief  summary  and  list  the  published  papers  which  describe 
the  completed  studies.  The  unpublished  results  are 
discussed  in  greater  detail. 

At  the  end  of  the  report  a  list  of  published  papers 
and  the  titles  of  graduate  theses  supported  under  this  contract 
is  given. 
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II.  WAVE  PROPAGATION  STUDIES 

In  this  area  the  efforts  have  been  concentrated  in 
studying  seismic  wave  propagation  in  heterogeneous  media. 

In  the  past  much  has  been  learned  about  the  earth's  interior 
and  the  seismic  source  from  the  interpretation  of  data  in 
terms  of  laterally  homogeneous  models  where  properties 
depended  only  on  depth.  Increased  data  on  the  lateral 
variations  of  properties  and  structure,  not  only  in  the 
crust  but  also  in  the  mantle,  have  made  it  necessary  to 
consider  seismic  surface  and  body  wave  propagation  in 
heterogeneous  media. 

By  the  nature  of  the  problem,  such  studies  can  best 
be  carried  out  by  numerical  methods.  Finite-difference 
techniques,  perturbation  methods,  and  statistical  techniques 
were  utilized  in  these  investigations.  Problems  of  interest 
included  surface  wave  propagation  in  media  with  non-parallel 
boundaries,  body  waves  incident  on  non-plane  boundaries, 
and  scattering  from  randomly  distributed  heterogeneities. 

For  the  propagation  of  body  and  surface  waves  in  media 
with  non-plane  boundaries  (such  as  sedimentary  basins  and 
the  crustal  structure  under  the  Large  Aperture  Seismic  Array 
(LASA)  m  Montana),  finite-difference  techniques  were  used. 
For  SH  waves  and  Love  waves  the  effects  of  heterogeneities 
were  found  to  be  significant.  Although  the  results  could 
have  been  predicted  intuitively  for  smooth  structures,  for 
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models  with  sharp  boundaries  or  heterogeneities  of  the  order 
of  a  wave  length  the  effects  become  very  complicated  because 
of  scattering ,  wave  conversion,  and  interference  phenomena. 

For  seismic  sources  in  regions  of  downgoing  slabs, 
both  ray  theory  and  wave  theory  formulations  have  been 
developed.  Very  strong  asymmetry  in  radiated  seismic  body 
waves  are  observed.  Strong  travel -time  anomalies,  amplitude 
variations  with  azimuth  and  distance  from  earthquakes  and 
explosions  such  as  Longshot  and  Milrow  can  be  explained  by 
the  lateral  variation  of  velocities  in  and  around  the  down¬ 
going  slab  under  island  arcs. 

The  splitting  of  the  eigenfrequencies  of  torsional 
free  oscillations  of  the  earth  due  to  lateral  heterogeneities 
is  studied  by  perturbation  theory.  The  perturbations 
are  expanded  in  spherical  harmonics  (0,<|>) ,  «ytyt  (0  ,<f>)  ] 

and  it  is  found  that  odd  terms  (s=odd)  produce  no  splitting 
o  first  order.  For  example,  there  would  be  no  splitting 
f  the  earth  had  a  continental  hemisphere  and  an  oceanic 
hemisphere.  Also,  heterogeneities  such  that  s>2£  will  not 
affect  the  eigenfrequencies  of  modes  of  order  A.  For  large 
1  this  means  that  free  oscillations  will  not  be  affected  by 
heterogeneities  with  scale  lengths  less  than  half  the  wave¬ 
length  of  the  free  oscillation.  The  amount  of  splitting  is 
computed  for  earth  models.  The  interpretation  of  high 
angular-order  modes  as  long-period  surface  waves  propagating 


along  great-circle  paths  are  analyzed  and  some  limits  are 
established  for  the  technique  of  phase  and  group  velocity 
averaging  along  great-circle  paths. 
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III.  SOURCE  MECHANISM  AND  FOCAL  DEPTH  STUDIES 

In  the  study  of  the  earthquake  source  mechanism  and 
focal  depth,  both  body  and  surface  waves,  with  special 
emphasis  on  the  latter,  have  been  used.  A  dislocation  model 
of  an  earthquake  source  was  used  for  determination  of  fault 
plane  geometry,  seismic  moment,  fault  length,  and  focal 
depth.  The  method  is  based  on  the  theory  of  normal  mode 
excitation  by  a  seismic  source  in  a  vertically  heterogeneous 
medium.  The  spectra  of  seismic  waves  recorded  at  distant 
stations  are  corrected  for  propagation  effects,  computed 
on  the  basis  of  structural  models.  Then  the  corrected  spectra 
are  interpreted  for  source  parameters. 

From  the  spectra  of  Rayleigh  and  Love  waves  at  a 
number  of  stations,  in  addition  to  the  dislocation  parameters 
and  fault  depth,  the  source  time  function  can  be  obtained. 

The  effects  of  these  parameters  on  the  surface  and  body 
wave  magnitudes  (Mg  and  m^) ,  and  hence  on  the  seismic 
discrimination  problem,  are  quite  significant.  Since  the 
source  dimensions  and  source  time  function  control  the  spectra, 
shapes,  they  also  determine  the  relative  amplitudes  of  body 
and  surface  waves  measured  at  1  Hz  and  0.05  Hz,  respectively. 
In  addition  to  above,  source  depth  is  an  important  parameter 
in  the  spectral  shape  of  surface  waves,  as  well  as  being 
a  discriminant  in  itself.  Thus  a  method  of  accurate 
determination  of  focal  depths  of  shallow  events  has  been 
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sought.  The  surface  waves  have  proved  to  be  most  sensitive 
to  focal  depth,  if  the  source  mechanism  and  source  time 
function  is  known,  or  kept  constant. 

Utilizing  Rayleigh  and  Love  waves,  source  characteristics 
of  both  underground  nuclear  explosions  and  many  earthquakes 
are  investigated.  These  are  discussed  in  great  detail  in 
the  publications  listed  below.  One  aspect  of  the  under¬ 
ground  explosion  source  function  which  complicates  the  above 
procedure  is  the  release  of  tectonic  strain  energy  by  the 
explosions.  Even  though  strong  Love  waves  are  generated  by 

some  explosions,  the  Mg  -  discriminant  is  not  seriously 
effected. 

Dislocation  models  were  applied  to  an  ensemble  of 
earthquakes  by  invoking  a  "similarity"  of  large  and  small 
earthquakes.  A  scaling  law  of  seismic  spectra,  derived  on 
the  basis  of  an  "w-squared"  model,  have  been  used  in  the 
discussion  of  various  problems  including  the  theoretical 
casis  for  the  Mg  -  m^  discriminant. 
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IV.  EARTH  STRUCTURE  AND  PROPERTIES 

The  elastic  and  anelastic  properties  of  the  earth  are 
both  radially  and  laterally  heterogeneous.  A  great  deal 
of  effort  has  gone  into  examining  the  detailed  variations 
of  these  properties,  with  particular  emphasis  on  understanding 
the  causes  of  the  variations  and  their  relationship  to  other 
geophysical  measurements  and  to  tectonic  patterns.  Dis¬ 
cussion  of  this  work  is  divided  below  into  three  sections: 

(1)  gross  earth  structure,  including  seismic  velocity  and 
density,  (2)  lateral  heterogeneities,  including  velocity, 
density,  and  Q,  and  (3)  the  important  effects  of  these 
findings,  particularly  those  on  Q,  on  seismic  wave  magnitudes 
and  the  discrimination  problem. 

1)  Earth  structure 

Various  techniques  have  been  used  to  infer  the  variation 
with  depth  of  seismic  velocity  and  density  in  the  earth. 

With  a  large  array,  the  slope  dt/dA  of  the  travel-time 
curve  can  be  measured  directly.  Using  LASA,  observations 
of  dt/dA  have  been  inverted  to  obtain  P- and  S-wave  velocities 
in  the  lower  mantle  (Chinnery  and  Toksfiz,  1967;  Toksflz 
ej:  al . ,  1967;  Fairborn,  1968,1969)  and  estimates  of  P-wave 
velocity  in  the  upper  mantle  (Wolfe,  1969).  An  important 
finding  was  the  discovery  of  several  sharp  increases  in 
velocity  with  depth  in  the  lower  mantle,  suggesting  that 
the  lower  mantle  is  not  homogeneous. 
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With  large ,  efficient  computers,  the  generation  of 
velocity-density  models  for  the  earth  that  satisfy  the 
entire  range  of  seismic  data,  as  well  as  mass  and  moment  of 
inertia,  is  a  simple  process.  The  problem  is  inherently 
non-unique,  however.  A  successful  attack  on  the  question 
of  non-uniqueness  is  the  Monte  Carlo  method,  in  which 
randomly  generated  earth  models  are  tested  against  obser¬ 
vations.  The  Monte  Carlo  technique,  which  offers  freedom 
from  bias  and  which  emphasizes  the  traits  common  to  successful 
models,  has  been  extensively  applied  to  the  problem  of 
earth  structure.  Highlights  of  the  conclusions  include 
constraints  on  the  composition  of  mantle  and  core  and 
insights  into  the  nature  and  thickness  of  lithosphere  and 

asthenosphere  (Press,  1968, 1969, 1970a, b;  Press  and  Kanamori, 
1970)  . 

2 )  Lateral  variations 


Both  the  elastic  and  anelastic  properties  of  the  crust 
and  upper  mantle  vary  with  tectonic  setting.  The  phase  and 
group  velocities  of  long-period  surface  waves  are  different 
for  oceans,  shields,  and  tectonically  active  regions. 
Inversion  of  the  dispersion  data  requires  that  significant 
differences  in  shear  velocity  among  the  various  regions 
persist  to  depths  of  at  least  400  km  (Toksflz  et  al. ,  1967). 

Even  more  regionally  variable  is  the  seismic  quality 
factor  Q,  particularly  within  the  low-Q,  low-velocity 
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asthenosphere.  These  lateral  variations  may  be  observed 
both  in  the  differential  attenuation  of  teleseismic  P-and 
S  waves  (Solomon  and  Toksflz,  1970)  and  in  the  attenuation 
of  surface  waves  of  period  sufficiently  long  so  that  a 
sizeable  fraction  of  the  wave  energy  propagates  in  the 
asthenosphere  (Tsai  and  Aki,  1969;  Solomon,  1971a). 

Because  the  asthenosphere,  at  least  in  tectonically 
active  areas,  is  almost  certainly  partially  melted, 
theoretical  considerations  suggest  that  Q  at  that  depth 
should  be  a  function  of  frequency.  This  is,  in  fact,  required 
by  the  data  for  western  North  America:  attenuation  and 
velocity  measurements  of  various  sorts,  in  all  spanning 
three  decades  in  frequency,  can  be  fit  only  by  postulating 
a  Q  that  depends  both  on  depth  and  frequency  (Solomon,  1971a, b) 
The  precise  dependence  of  Q  on  frequency  xs  governed  by 
properties  of  the  partial  melt  such  as  melt  viscosity  and 
volume  concentration  of  melt.  All  of  the  lateral  variation 
of  Q  in  the  asthenosphere  may  be  explained  by  relatively 
modest  lateral  changes  in  temperature  (Solomon,  1971b). 

Among  the  most  dramatic  lateral  variations  in  such 
quantities  as  velocity  and  Q  occur  at  island  arcs,  where  the 
relatively  cool  lithosphere  is  being  undsrthrust  into  the 
mantle.  Seismic  velocities  calculated  from  realistic 
temperature  models  for  downgoing  lithospheric  slabs  (Minear 
and  Toksflz,  1970a, b;  Toksflz  et  al. ,  1971)  predict  that  the 
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slab  will  strongly  affect  the  travel-time  and  propagation 
direction  of  body  waves  from  nearby  sources.  This  has 
important  implications  for  determination  of  location  and 
source  mechanism  for  seismic  events  in  such  regions  (Toksflz 

et  al.  1971)  . 

Spherical  harmonic  expansion  and  correlation  of  such 
globally  distributed  geophysical  data  as  seismic  travel-time 
anomalies,  heat  flow,  surface  topography  and  gravity 
demonstrate  both  the  interrelationships  of  these  various 
quantities  and  the  nature  of  the  global  variations  (Toksflz 
and  Arkani-Hamed,  1967;  Arkani-Hamed  and  Toksflz,  1968; 

Toksflz  et  al. ,  1969;  Arkani-Hamed,  1969,1970).  A  major 
conclusion  of  such  studies  is  that  lateral  variation  of  such 
properties  as  density  extends  well  into  the  lower  mantle. 

3)  Some  effects  on  magnitudes  and  discrimination 

The  nature  of  heterogeneities  in  the  earth  must  be 
well  understood  before  seismic  wave  amplitudes  may  be 
corrected  for  propagation  effects  to  estimate  source  parameters. 
Foremost  among  the  relevant  earth  properties  is  Q,  which 
controls  both  the  absolute  amplitudes  and  spectral  shape  of 
seismic  waves. 

Lateral  variation  of  Q  can  appreciably  shift  the 
MS  ~  mb  relati°nship.  In  particular,  in  tectonically  active 
areas  such  as  parts  of  western  North  America,  mb  is  generally 
lowered  by  several  tenths  of  a  magnitude  unit,  relative  to 
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events  of  similar  Mg  in  more  stable  regions,  due  to 
abnormally  high  attenuation  in  the  upper  mantle  of  the 
tectonic  region  (Solomon  et  al .  ,  1970;  Ward  and  Toksfiz,  1971). 
The  reason  is  simple.  Because  surface  waves  in  the  period 
range  used  to  define  Mg  propagate  principally  in  the  high-Q 
lithosphere,  and  because  scattering  at  ocean-continent 
boundaries  is  generally  small  at  these  periods,  determination 
of  Mg  is  relatively  insensitive  to  path.  For  body  waves 
that  have  propagated  through  the  low-Q  asthenosphere  in  which 
Q  is  strongly  variable  from  one  region  to  another,  m^  is 
strongly  path  dependent. 

In  any  amplitude  measurements  that  cover  a  broad 
frequency  band,  some  account  needs  to  be  made  of  the 
regional  variations  and  frequency  dependence  of  0  in  order 
to  extract  the  maximum  possible  information  about  the  source!. 
In  particular,  determination  of  source  depth  from  surface 
wave  spectral  shape  (see  section  III)  and  of  source 
dimensions  and/or  time  function  from  body  wave  spectra  must 
be  based  on  an  accurate  knowledge  of  Q  in  the  earth. 
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V.  SEISMIC  ARRAY  STUDIES 

With  the  installation  of  the  Large  Aperture  Seismic 
Array  (LASA) ,  a  new  capability  was  added  to  seismological 
studies.  Availability  of  seismic  data  in  digital  form  from 
a  series  of  sensors  has  opened  new  areas  of  study.  Under 
this  contract  we  have  been  involved  in  three  separate 
investigations:  (1)  Measurement  of  apparent  slowness  of 

body  waves  (dt/d  A  )  and  the  interpretation  of  these  in 
terms  of  detailed  earth  structure,  (2)  Frequency-wave 
number  study  of  the  nature  of  seismic  noise  and  microseisms, 
and  (3)  Evaluation  of  the  automatic  event  detection 
capability  of  arrays. 

1)  dt/dA  measurements 

The  dt/dA  measurements  for  seismic  P  and  S  waves  have 
greatly  enhanced  our  ability  to  determine  the  finer  features 
of  the  earth's  mantle.  These  results  were  discussed  in 
greater  detail  in  the  previous  section. 

2)  Frequency-wave  number  spectra  and  nature  of  microseisms 

Utilizing  data  from  LASA,  microseisms  have  been  studied 
by  computing  the  frequency-wave  number  spectra.  Wave  types, 
modes  of  propagation  and  source  regions  have  been  identified. 
At  LASA  shoirt-period  (0.5  -  2.0  sec)  microseisms  consist 
primarily  of  coherent  body  waves,  generated  at  the  oceanic 
storm  centers.  At  intermediate  periods  (2-5  sec)  higher¬ 
mode  surface  waves  dominate.  At  still  longer  periods,  the 
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microseisms  are  made  of  fundamental  mode  Rayleigh  waves  and 
some  Love  waves . 

3)  Automatic  event  detection 

Since  April  1969  detection  of  seismic  signals  has  been 
done  automatically  and  cn-line  by  the  deployment  of  300 
pre-steered  beams.  Each  tine  the  signal  to  noise  ratio 
goes  above  the  10  db  threshold  in  any  beam  the  event 
detector  records  a  detection  and  groups  it  together  with 
any  other  beam  detections  occurring  at  the  same  time.  The 
output  of  the  detector  is  a  list  of  detection  groupings  with 
their  times,  intensity,  and  beam  number  of  maximum  intensity. 
There  are  about  160  detection  groups  a  day  of  which  only 
25  are  genuine  signals.  This  study  was  undertaken  to  find 
how  one  can  separate  the  signals  from  the  false  alarms  and 
how  one  can  classify  the  signals  into  first  arrivals,  P  or 
PKP  and  the  later  phases  PcP,  ScP,  etc.  using  some  auto¬ 
matic  scheme. 

If  processing  of  earthquake  data  could  be  done  auto¬ 
matically,  then  earthquake  reporting  would  become  more  uniform 
and  standardized  with  time.  Furthermore,  more  phases  would 
be  picked  up  since  the  weak  phases  would  not  be  ignored  if 
they  happened  to  fit  travel-time  tables  and  come  from  the 
same  azimuth  as  the  earthquake.  If  seismic  data  is 
available  only  from  one  localized  array  then  information  from 
later  phases  is  very  useful  in  improving  estimates  of  the 
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earthquake  location.  In  some  cases  the  detection  of  a  later 
phase  would  confirm  the  interpretation  of  a  first  arrival. 
Later  ohases  may  contain  relevant  information  about  the 
focal  mechanism. 

In  the  first  figure  the  percent  detection  by  LASA  of 
Coast  and  Geodetic  Survey  (CGS)  reported  events  of 
magnitude  m^  and  greater  is  plotted  for  two  distance  ranges. 
The  earthquakes  missed  when  the  LASA  system  was  not  oper¬ 
ating  wore  not  taken  into  account.  The  number  of  events 
missed  in  the  distance  range  30  to  90  degrees  appears  to  be 
independent  of  magnitude.  Since  from  experience  it  is  known 
to  be  highly  unlikely  for  LASA  to  miss  an  event  of  magnitude 
4.5  or  larger  in  this  distance  range,  it  can  be  concluded 
that  LASA  has  the  capability  to  detect  all  events  that  CGS 
lists  in  the  30  to  90  degrees  distance  range.  Furthermore, 
there  are  many  more  events  that  LASA  reports  but  are  unlisted 
in  the  CGS  bulletins.  At  larger  distances  near  and  beyond 
the  shadow  zone,  LASA  misses  not  more  than  40%  of  the 
events  reported  by  CGS.  A  magnitude  dependence  can  be  seen. 
Thus  it  follows  that  though  a  single  array  does  not  excel 
in  location  capabilities  it  certainly  has  exceptionally 
good  detection  capabilities. 

In  developing  the  automatic  phase  identifier  the 
guiding  principle  was  to  use  the  same  ideas  and  methods  that 
an  analyst  would  employ.  Given  any  pair  of  detections  the 
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parameters  of  the  detections,  namely  velocity,  azimuth, 
signal  intensity,  and  time  interval  between  detections  must 
satisfy  a  certain  relationship  (with  some  statistical  error) 
in  order  to  be  classified  as  a  certain  pair  of  phases 
related  to  the  same  earthquake. 

Phase  Identifier  -  The  input  to  the  phase  identifier  is  a 
list  of  detection  times  in  chronological  order  with  their 
signal  strength  and  beam  number  of  maximum  intensity.  The 
output  is  the  interpretation  of  these  detections,  p,  PKP, 
noise  or  later  phase.  The  phase  identifier  will  recognize 
the  following  later  phases:  PcP,  ScP,  PP,  SKS,  PKKP,  and 
P'P'.  In  addition  it  will  find  the  first  arrival  associated 
with  the  later  phase.  The  phase  identifier  cannot  pick  out 
a  later  phase  in  the  detection  log  it  the  first  arrival 
associated  with  the  later  phase  was  missed.  Thus  it  tends 
to  miss  some  later  phases  coming  from  events  around  110 
degrees  because  the  first  arrival  was  undetected. 

Due  to  the  fact  that  the  input  information  is  inexact, 
the  statistical  approach  was  used  in  classifying  a  signal. 

The  design  criterion  was  to  minimize  the  misclassif ication 
error  hence  the  optimum  scheme  is  the  maximum  likelihood 
test.  Using  a  set  of  detections  in  which  the  classification 
of  the  detections  was  known  from  outside  sources  such  as  the 
CGS  catalog  and  Seismic  Array  Analysis  Bulletins  (SAAC) ,  an 
estimate  was  made  of  the  likelihood  function  given  a  set  of 


23. 


observations  and  a  specific  interpretation. 

Ai  =  Pr  (NBMl , MSTAl , AT ,  NBM2 ,MSTA2 1  Hi )  (1) 

is  the  probability  of  observing  the  first  arrival  in 
beam  NBMl  with  intensity  MSTAl  followed  by  a  later  arrival 
AT  seconds  later  in  NBM2  and  with  intensity  MSTA2  assuming 
that  the  detection  pair  are  identified  by  hypothesis 
(e.g.  PKP-PKKP) .  The  likelihood  function  was  approximated 
and  stored  by  the  phase  identifier  so  that  gi /en  any 
detection  pair  with  their  parameters  and  given  an  inter¬ 
pretation,  the  likelihood  of  observing  these  parameters 
for  that  interpretation  could  be  estimated.  The  likelihood 
function  is  computed  for  all  possible  detection  pairs 
separated  by  less  than  30  minutes  and  for  all  possible  inter¬ 
pretations,  and  then  the  interpretation  with  the  largest 
likelihood  is  chosen.  The  bulk  of  this  section  will  deal 
with  the  problem  of  how  to  estimate  A^. 

It  is  not  practical  to  work  with  the  parameters  NBMl, 
MSTAl,  AT,  NBM2 ,  MSTA2  since  they  are  generally  not  independent 
of  each  other.  It  is  very  desirable  if  A^  could  be  approx¬ 
imated  by  a  set  of  separable  parameters,  i.e. 


=  Pr(Sl|Hi)Pr(S2|  Hi)Pr(S3|Hi) 


(2) 
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where  PrfSllH^)  is  the  probability  of  observing  parameter 
Si  given  the  hypothesis  lb.  For  the  set  of 'new  parameters 

the  following  were  adopted: 

51  -  DIS(NBMl)  -  DIS(NBM2jH.) 

52  =  DIS(NBMl)  -  CIS ( At|  lb  ) 

53  =  AZ(NBMl)  -  AZ(NBM2|lb) 

MSTAl  =  MSTAl 
AMSTA  =  MSTAl  -  MSTA2 

where  DIS(NBM2|H^)  is  the  distance  associated  with  NBM2 
given  it  is  some  phase  lb,  and  AZ  is  the  azimuth  associated 
with  the  respective  beam  and  hypothesis.  DIS(AT|lb)  is  the 
distance  inferred  from  travel-time  tables.  (Depth  corrections 
were  unnecessary;  see  Figures  2  and  3.)  These  parameters 
should  be  approximately  independent  of  each  other. 

These  parameters  have  another  important  property.  if 
the  specific  hypothesis  II.  is  correct  then  the  two  detections 
are  associated  with  the  same  event  and  SI,  S2,  and  S3  should 
cluster  around  zero.  On  the  other  hand  if  the  interpretation 
is  wrong  then  SI,  S2,  and  S3  can  take  any  value. 

In  actual  cases  SI,  S2  and  S3  can  have  a  considerable 
departure  from  zero  under  a  correct  interpretation.  Figure 
4  plots  CGS  epicenters  of  specific  phases  comir.n  in  specific 
beams.  The  scatter  can  be  attributed  to  incc  eroncies 
among  the  waves  from  the  different  subarray  and  th.  fact  that 
with  a  finite  number  of  beams  it  is  unlikely  for  the  chase 


25. 


to  fall  directly  in  a  beam.  The  distances  and  azimuths 
associated  with  the  transformation  were  determined  using  a 
training  set  of  detections,  a  set  of  2000  detections  which 
vere  identified  from  CGS  and  SAAC  bulletins.  There  were 
many  beams  and  phases  for  which  no  detections  were  found 
in  the  training  set.  For  instance  there  were  only  30 
beams  in  which  P'P*  detections  were  observed.  If  for  a 
specific  phase  no  detection  in  the  training  set  was  observed 
in  the  beam  then  that  phase  was  categorically  rejected  in 
the  present  implementation  of  the  phase  identifier.  Thus 
the  phase  identifier  would  tend  not  to  classify  phases 
coming  from  aseismic  regions.  (This  restriction  could  be 
removed  for  strong  signals.) 

In  the  upper  half  of  Figure  5  the  marginal  distribution 
of  SI  is  shown  for  correctly  identified  phases  shown  in 
the  lower  half.  (A  synthetic  detection  log  described  in  the 
next  section  was  used  to  get  the  lower  distributions.)  The 
lower  distributions  are  not  uniform  since  there  was  a 
certain  amount  of  screening  before  considering  an  interpre¬ 
tation.  For  example  the  phase  identifier  would  not  consider  the 
P-P’P'  interpretation  if  AT  was  less  than  1400  seconds. 

These  distributions  summarize  the  usefulness  of  the 
parameters  SI,  S2,  and  S3.  SI  and  S2  are  necessary  to 
distinguish  the  later  phases  since  they  are  the  only  parameters 
which  depend  strongly  on  the  interpretation  (S3  only  separates 
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PcP,  ScP,  SKP,  PP  interpretations  from  PKKP  and  P'P’). 

S3  shows  the  best  contrast  under  the  two  hypotheses  and  is 
very  useful  in  deciding  that  the  two  signals  came  from  the 

same  earthquake. 

The  covariance  matrix  for  SI,  S2  and  S3  is  shown  in 
Table  1  for  correct  and  incorrect  identifications.  The 
cross  terms  were  small  enough  to  neglect. 

The  remaining  parameters,  MSTAl  and  MSTA  are  principally 
used  to  distinguish  a  genuine  signal  from  a  false  alarm. 

The  MSTA  is  the  maximum  short  term  average,  averaged  over 
1.8  seconds  of  filtered,  rectified  aliened  beam  data  sampled 
10  times  a  second.  It  is  measured  in  so  called  quantum 
units.  In  Figure  6  MSTA  is  plotted  versus  magnitude  of 
shallow  events  at  a  distance  between  30  and  60  degrees.  (The 
attenuation  factor  is  fairly  constant  in  this  range.) 

There  appears  to  be  a  better  correlation  suggesting  that 
MSTA  is  probably  more  logarithmically  related  to  the  signal 
amplitude;  however,  the  scatter  is  so  large  that  one  should 
not  attempt  to  estimate  the  magnitude  from  MSTA. 

In  Figure  7  the  cumulative  frequency  of  detections 
versus  MSTA  is  plotted  for  all  detections  in  the  time 
interval  April  13,  1969  to  October  1,  1969.  in  the  range 
MSTA  between  400  and  1000  quantum  units,  the  log  frequency 
distribution  is  linear.  Above  1000  the  signals  in  the 
different  subarrays  are  clipped. 


For  small  MSTA  the  number 
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of  detections  increases  abnormally.  It  was  observed  that  the 
strong  signals  seem  to  be  confined  to  a  limited  set  of 
beams  corresponding  to  seismic  regions  while  the  weak  signals 
occur  uniformly  over  both  seismic  and  eiseismic  regions.  In 
Figure  8  cumulative  distribution  of  detections  is  plotted 
versus  MSTA  for  two  sets  of  36  beams,  the  first  set  has  61 
percent  of  the  detections  with  MSTA  >  250  while  the  second 
set  has  less  than  one  percent  of  the  detections  with  MSTA  >  250. 
The  distributions  are  very  different.  It  seems  reasonable 
to  conclude  from  these  results  that  the  linear  part  of  the 
log  frequency  distribution  in  Figure  7  is  due  to  real 
earthquakes  and  probably  reflects  Gutenberg-Richter ' s 
frequency-magnitude  distribution;  the  nonlinear  portion  of 
the  distribution  probably  is  not  due  to  real  earthquakes 
but  are  false  alarms  occurring  when  the  0.9  to  1.4  hz  noise 
component  becomes  suddenly  larger  in  some  specific  beam. 

It  appears  that  the  frequency  of  false  alarms  is  inversely 
proportional  to  the  fourth  power  of  MSTA.  In  Figure  9  the 
cumulative  false  alarm  percent  of  the  total  number  of 
detections  is  plotted  versus  detections  with  MSTA  above  a 
variable  threshold. 

So  far  the  MSTA  distributions  have  been  determined  for 
genuine  signals  and  noise.  The  next  parameter  MSTA2  would 
depend  on  MSTA1  if  both  detections  are  related  to  the  same 
earthquake.  Assuming  that  MSTA  is  logarithmically  related 
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to  signal  strength  then  AMSTA=MSTA1-MSTA2  would  be  related 
to  the  attenuation  factor.  Prom  earlier  studies  it  has  been 
observed  that  the  attenuation  factor  is  very  variable  even 
for  the  same  phase  coming  from  the  same  region.  The  only 
consistent  observation  is  that  later  phases  are  rarely 
Stronger  than  the  first  arrival  except  in  the  shadow  tone. 

The  distribution  of  AMSTA  was  estimated  and  approximated 
by  two  exponentials. 

At  this  point  there  is  sufficient  information  available 
to  estimate  Aj.  using  the  transformed  set  of  parameters  and 
equation  (2,.  Table  2  lists  the  actual  distributions  used 
to  approximate  the  probability  density  functions  of  si,  S2, 
S3,  MSTA1  and  AMSTA.  Time  was  saved  by  computing  log  A. 
since  this  is  just  a  monotonic  increasing  transformation. 

In  summary,  the  phase  identifier  works  as  follows 
All  detections  that  have  occurred  in  the  last  half  hour  and 
that  have  triggered  a  beam  from  which  P  or  PKP  phases  were 
observed  in  the  training  set  are  stored  in  a  memory  buffer. 

Any  new  detection  is  then  tested  with  each  of  the  earlier  ' 
detections  in  the  buffer.  For  a  specific  later  phase  to  be 
considered  AT  of  the  detection  pair  must  fall  i„  a  specifio 
time  window  of  that  interpretation.  Next  the  phase 
identifier  would  attempt  to  transform  to  new  variables  SI, 

If  the  later  phase  was  never  observed  in  NBM2, 
then  that  interpretation  would  be  rejected  at  this  point. 
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LogfiA  would  be  computed  by  equation  (3) 

l°geA  “  lo9eAi  "  lo9eAf  (3) 

where  Af  is  computed  assuming  the  interpretation  is  false. 

If  logeA  is  negative  it  means  that  it  is  more  probable  that 
the  interpretation  was  incorrect  and  hence  it  is  immediately 
rejected.  On  the  other  hand  if  it  is  positive  it  is  now 
ascertained  that  the  new  detection  is  a  later  phase.  The 
final  interpretation  for  the  new  phase  is  the  one  which 
maximizes  logeA  .  If  all  possible  interpretations  were 
rejected  then  the  phase  identifier  concludes  the  new  de¬ 
tection  is  not  any  of  the  later  phases  it  was  designed  to 
identify.  Hence  the  new  detection  is  either  a  P,  PKP,  or 
noise,  if  no  first  arrivals  were  observed  in  NBM2  then  the 
detection  is  taken  to  be  noise,  otherwise  it  is  tentatively 
assumed  to  be  a  first  arrival  and  stored  in  the  buffer. 

Unless  the  MSTA  happened  to  be  large,  one  could  never  know 
for  sure  whether  the  new  arrival  was  noise  or  a  first  arrival. 
The  only  exception  is  when  a  later  arrival  is  discovered  to 
correspond  to  the  first  arrival.  As  will  be  described  in 
the  next  section  this  implementation  of  phase  identifier 
appears  to  work  very  well. 

Performance  -  The  phase  identifier  described  in  the  last 
section  was  run  through  67  days  of  actual  data.  On  the 
average  it  picked  out  9  later  phases  a  day.  Of  these  9 
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later  phases  only  three  of  them  could  be  confirmed  by  using 
outside  bulletins.  About  half  of  these  phases  are  so  small 
that  an  analyst  would  probably  not  notice  them.  There  was 
no  way  by  which  we  could  confirm  the  remaining  6  phases 
since  corresponding  earthquakes  were  not  found  in  the  SAAC 
or  CGS  catalogs.  However,  it  shall  be  shown  that  it  is 

very  probable  that  3  of  these  unconfirmed  later  phases  were 
correctly  identified. 

The  phase  identifier  made  few  errors  in  misclassifying 
a  later  phase.  Later  phase  identifications  in  the  SAAC 
bulletins  were  compared  with  those  put  out  by  the  phase 
identifier.  Later  phases  in  the  SAAC  bulletin  that  did  not 
cause  a  detection  were  not  included  in  this  analysis.  The 
results  are  shown  by  the  confusion  matrix  in  Table  3.  The 
row  is  the  SAAC  identification.  Along  any  row  is  the 
distribution  of  phase  identifier  decisions.  The  sum  of  the 
otf  diagonal  terms  (misclassif ications)  is  only  10%  of  all 
the  SAAC  identifications.  The  missed  phases  were  due  generally 
to  no  first  arrival  in  the  detection  log  either  because  the 
signal  was  too  weak  or  the  LASA  detection  system  was  down. 

In  some  cases  later  phases  were  missed  because  the  phase 
identifier  was  not  expecting  to  see  the  phase  in  that 
particular  beam.  The  later  phases  that  were  not  missed  but 
misidentified  were  due  to  the  confusion  resulting  when  two 
earthquakes  occurred  around  the  same  time  and  along  the  same 
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azimuth. 

It  is  still  necessary  to  estimate  the  probability  of 
the  phase  identifier  picking  a  later  phase  which  is  actually 
noise  or  a  first  arrival.  Since  it  is  impossible  to  know 
from  outside  sources  that  a  specific  detection  is  noise 
another  approach  had  to  be  taken.  The  identifying  process 
is  so  complicated  that  it  is  very  difficult  to  estimate 
the  number  of  errors  using  a  theoretical  model.  The  most 
reliable  way  and  also  the  easiest  way  is  by  simulation  using 
a  synthetic  detection  log. 

Using  a  random  number  generator  a  new  detection  log 
was  generated  with  all  the  statistical  properties  of  the 
original  detection  log  except  one.  Thare  were  no  later 
phases  in  the  synthetic  detection  log,  i.e.  all  the  de¬ 
tections  were  completely  independent  of  each  other.  The 
method  used  in  synthesizing  this  log  is  described  in  Appendix 
A.  When  this  log  was  run  through  the  phase  identifier  less 
than  two  later  phases  were  picked  out  a  day  on  the  average. 
Thus  it  could  be  concluded  that  80«  of  the  later  phases 
found  by  the  phase  identifier  are  true  later  phases. 

At  the  present  the  automated  phase  identification  and 
event  verification  is  being  extended  to  two  arrays  (LASA 
and  NORSAR)  operating  simultaneously.  These  provide  better 
geographical  coverage,  greater  reliability  and  eventually 
better  automatic  locations  of  events. 


32. 


Publications 


Chinnery,  M.A.  and  M.N.  Toksfiz,  P-wave  velocities  in  the  mantle 
below  700  km,  Bull.  Seism.  Soc.  Am.,  5.7,  199-226,  1967. 

Lacoss,  R.T.,  E.J.  Kelly  and  M.N.  Toksdz,  Estimation  of  seis¬ 
mic  noise  structure  using  arrays,  Geophysics,  21-38, 

1969. 

Shlien,  S.  and  M.N.  Toksflz,  A  clustering  model  for  earthquake 
occurrences,  Bull.  Seism.  Soc.  Am.,  60,  1765-1787,  1970. 

Shlien,  S. ,  and  M.N.  Toksfiz,  Frequency-magnitude  statistic 
of  earthquake  occurrences,  Earthquake  Notes,  4^,  5-18, 

1970. 

Toksfiz,  M.N.,  Crustal  effects  on  long-period  chirp  filters, 
Proceedings  of  the  Woods  Hole  Conference  on  Seismic 
Discrimination,  Vol.  I.  20-23  July  1970. 

Toksfiz ,  M.N. ,  M.A.  Chinnery  and  D.L.  Anderson,  Inhomogeneities 
in  the  earth's  mantle,  Geophys.  J.  R.  astr.  Soc.,  13, 

31-59,  1967. 

Toksfiz ,  M.N.  and  R.T.  Lacoss,  Mode  structure  and  source  of 
microseisms,  Science,  159,  872-873,  1968. 

Toksfiz ,  M.N.  and  R.A.  Wiggins,  Seismic  arrays  and  the 

structure  of  the  earth's  interior,  NEREM  Record,  170-171, 
1969. 

Wolfe,  J.C.,  dt/dA  measurement  of  short  period  P  waves  in  the 
upper  mantle  for  the  Western  North  America  using  LASA, 

S.M.  thesis,  M.I.T.,  1970. 


TABLE  1 


Covariance  Matrix 


SI 

S2 

S3 

156 

34 

18 

34 

88 

13 

18 

13 

236 

2000 

103 

210 

103 

499 

-57 

210 

-57 

10650 

TABLE  2 . 


N(a,b)  Normal  Distribution 
Mean  a,  Variance  b 


True 

False 

SI 

N (0 , 111) 

N  ( 10 . 5 , 16  70 ) 

S2 

N (0 , 72) 

N  (1. 4 , 728) 

S3 

N (0,202) 

1/360 

MSTA1 

1  -MSTA1/420 

420  e 

375,000/ (MSTA1 ) 4 

AMSTA 

.88  -AMSTA/580 

580 

AMSTA  >  0 

.12  AMSTA/80 

88  e 

AMSTA  <  0 

375,000 

TMSTA1  4-  AMSTA)  T 

A. A 
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TABLE  3. 


Confusion  Matrix 


Phase  Identifier 

PCP 

SCP 

PP 

SKP 

PKKP 

p.p. 

MISSED 

PCP 

39 

1 

0 

0 

0 

0 

3 

SCP 

0 

11 

1 

1 

1 

0 

0 

pp 

0 

0 

16 

0 

0 

0 

4 

SKP 

0 

0 

0 

11 

0 

0 

2 

PKKP 

0 

0 

0 

0 

38 

0 

0 

pipi 

0 

0 

0 

0 

0 

6 

0 
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FIGURE  CAPTIONS 

5-1.  Percent  detection  of  CGS  events  magnitude  m}  and 
greater.  Right  events  between  31  and  90  degrees 
from  LASA;  left  events  between  90  and  150  degrees 
from  LAS A. 

5-2.  Time  interval  between  later  arrival  and  first  arrival 
in  seconds  versus  CGS  distance  from  LASA  in  degrees 
for  PcP,  ScP,  PP-P,  PP-PKP ,  SKP  phases. 

5-3.  Time  interval  between  later  arrival  and  first  arrival 
in  seconds  versus  CGS  distance  from  LASA  in  degrees 
for  P'P',  PKKP-P ,  PKKP-PKP  phases. 

5-4.  Epicenter  location  with  respect  to  LASA  for  phases 
arriving  in  specified  beams. 

5-5.  Frequency  distribution  of  SI  in  degrees.  Top  are 

phases  picked  out  by  the  phase  identifier  (657  total 
samples) ;  bottom  phases  from  synthetic  detection  log 
(733  samples). 

5-6.  MSTA  of  detections  of  P  phase  in  quantum  units  versus 
CGS  magnitude.  Earthquakes  were  shallow  and  confined 
to  30  to  60  degrees  from  LASA. 

5-7.  Cumulative  number  of  detections  versus  MSTA  in  quantum 
units  for  141  days  from  13  May  to  1  October  1969. 

19253  detections  with  MSTA  >  50. 
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5-8.  cumulative  number  of  detections  versus  MSTA  in  quantum 
units.  Upper  distribution  determined  from  36  beams 
containing  61%  of  detections  with  MSTA  >  250  quantum 
units.  Lower  distribution  from  36  beams  containing 

less  than  0.65  of  detections  with  MSTA  >  250  over  the 
same  time  interval. 

3-9.  Estimated  false  alarm  rate. 
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/appendix  A 

In  this  section  we  describe  how  the  synthetic  detection 
log  was  generated.  The  synthetic  detection  log  had  to 
have  the  following  properties.  (1)  Detections  occur  at 
completely  random  times,  (2)  About  140  detections  a  day 
have  a  signal  strength  distributed  as  MSTA-4  and  occur 
uniformly  in  all  300  beams,  (3)  About  25  detections  a  day 
are  confined  to  seismic  beams  and  have  an  exponential  MSTA 
distribution. 

The  random  number  generator  produced  uncorrelated 

numbers  uniformly  distributed  between  plus  and  minus 
17 

2  =  131072  which  were  transformed  to  the  interval  0  to  1. 

To  satisfy  the  first  condition  it  was  assumed  that  detection 
occurrence  times  could  be  modeled  by  the  Poisson  process. 

Then  the  time  intervals  between  adjacent  detections  are 
independently  and  exponentially  distributed.  The  expo¬ 
nentially  distributed  time  intervals  t  in  seconds  were 
generated  by  the  following  equation  using  the  uniformly 
distributed  random  numbers  x,  between  0  and  1 

1  =  -  k  loge(1_x)  (1A) 

k  is  the  number  of  detections  in  a  day  divided  by  the  number 
of  seconds  in  a  day  (86,400). 

For  the  second  and  third  conditions  another  random 
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number  was  picked.  If  the  number  was  below  a  certain  threshold 
then  the  detection  was  taken  to  be  a  seismic  signal ,  other¬ 
wise  it  was  noise.  The  threshold  was  chosen  so  that  there 
were  about  25  seismic  signals  a  day.  If  the  seismic  signal 
was  noise  then 

.4 


Pr  (MSTA| noise)  =  A/MSTA^ 


MSTA  >50 


HA) 


where  the  normalization  constant  A  =  3*504  =  375000.  Another 
x  was  computed  and  MSTA  was  determined  by 

1/3 


MSTA  =  (A/3x) 


'  3A) 


The  beam  number  was  determined  so  that  it  was  uniformly 
distributed  between  1  and  300. 

If  the  signal  was  seismic  then  MSTA  probability  density 
function  is  given  by 

Pr(MSTA| seismic)  =  (1/420 ) *exp (-MSTA/420)  (4A) 


Hence  MSTA  was  determined  by 


MSTA  =  -4201ogp  (1-x) 


(5A) 


The  beam  number  was  determined  so  that  it  followed  the 

seismic  distribution  in  Table  1. 

It  was  found  that  the  number  of  later  phases  picked  out 

by  the  phase  identifier  was  fairly  insensitive  to  k  but 
depended  on  the  seismic  beam  distribution. 
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VI .  INSTRUMENTATION 

Long-period  tiltmeters  were  installed  at  two  stations  - 
Harvard,  Mass,  and  Eilat,  Israel  -  under  this  contract. 

These  tiltmecers  are  redesigned  versions  of  the  Benioff 
mercury  pendulum.  The  mechanical  stability  is  improved 
by  manufacturing  tanks  of  glass.  Bridge  electronics 
operates  on  a  modified  concept  with  improved  stability  and 
S/N  ratio.  An  inert,  stable  transformer  oil-fills  the  gap 
between  the  mercury  pool  surfaces  and  capacitor  plates. 

The  base  length  of  each  component  of  the  tiltmeter 
at  Harvard  is  19  ft,  and  the  natural  period  is  70  sec. 

Linear  tilt  range  is  2  x  10  radians.  Resolution  is  1  x  10~ 
radians.  At  the  Israel  installation  the  base  length  is 
increased,  resulting  in  greater  sensitivity  and  longer 
natural  period  (Tq  =  ,.02  sec)  .  The  resolution  is  1  x  lO-11 
radians.  At  both  installations  the  instruments  record  in 
three  pass-bands:  (1)  Surface  wave  band,  (2)  Eigen- 
frequency  band,  and  (3)  D.C.  output. 

Because  of  greater  sensitivity,  a  seismically  quieter 
site  (low  level  of  microseisms)  and  better  temperature 
stability  at  the  vault,  the  performance  of  the  tiltmeter 
installed  at  Eilat  is  far  superior  to  that  at  Harvard  station. 

As  of  January  1971,  about  three  months  of  records  were 
available  to  assess  the  performance  of  the  tiltmeters 
installed  in  Eilat.  This  is  not  yet  a  sufficient  suite  of 
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records  to  achieve  the  research  goals  referred  to  earlier, 
but  it  does  provide  preliminary  indications  on  the  per- 
formance  of  the  instruments. 

An  examination  of  the  records  indicates  that  in  the 
surface -wave  band  the  instruments  are  operating  near  the 
maximum  sensitivity  permitted  by  the  ground  noise  character¬ 
istic  of  the  winter  months,  widening  the  gap  of  the 
transducer  capacitance  will  lead  to  some  improvement  in 
performance  in  the  winter  months ,  and  to  a  marked  improvement 
in  the  summer  and  on  quiet  winter  days.  The  surface -wave 
channel  shouid  be  markedly  improved  by  this  modification. 

The  records  differ  from  conventional  seismographs  in 
showing  mantle  surface  waves  from  relatively  small  events. 

We  have  not  been  able  to  identify  most  of  these  events 
because  their  magnitudes  are  too  small  to  be  reported  through 
usual  channeis.  m  further  studies,  we  will  track  those 
tremors  down  in  order  to  study  the  efficiency  of  long  wave 
excitation  in  the  magnitude  range  less  than  5. 

Some  sample  records  are  shown  in  Fig.  6-1.  ln  examining 
new  records,  the  reader  should  remember  that  the  response 
at  20  seconds  is  down  by  a  factor  of  about  5  compared  to 
the  response  at  70  seconds,  so  that  performance  of  the 
surface -wave  channei  should  be  judged  by  the  detection  of 
waves  with  periods  in  the  range  50-100  seconds,  other 
seismographs  are  better  suited  for  recording  surface  waves 
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with  periods  under  30  seconds,  and  the  eigenmode  channel 

will  be  more  useful  for  waves  with  periods  exceeding  200 
seconds . 

The  records  in  Fig.  6-1  are  typical  of  a  large  number 
of  events  detected  in  the  three-month  operating  period 
showing  crustal-mantle  surface  waves  containing  energy 
near  100  seconds  from  earthquakes  with  magnitudes  near  5. 
The  background  noise  shown  in  Fig.  6-1  is  typical, 

indicating  that  the  gain  could  have  been  increased  by  two 
or  more  times  about  50%  of  the  time. 
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Record  1.  Prince  Edward  Island.  TV!  =5.1,  A  =  75  .  Showing  excitation 
80- second  surface  waves  on  E-W  component. 


Record  2.  California,  M«5.4.  A  =  110°.  Showing  excitation  of 
100-second  surface  waves  on  N-S  component. 


Record  3. 


Jon  Ma|cn  Island,  M  =  5.  1,  A=  47°.  Showing  excitation 
of  80  -second  surface  waves  on  K-\V  component. 


Fig.  6-1  Surface  wave  records 


*»* 
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